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ABSTRACT 

I present a general classification of self-similar solutions to the equations of 
gravitational hydrodynamics that contain many previous results as special cases. For 
cold flows with spherical symmetry, the solution space can be classified into several 
regions of behavior similar to the Bondi solutions for steady flow. A full description of 
these solution is possible, which serves as the asymptotic limit for the general problem. 
By applying a shock jump condition, exact general solutions can be constructed. The 
isothermal case allows an extra exact integral, and can be asymptotically analyzed in 
the presence of finite pressure. These solutions serve as analytic models for problems 
such as spherical accretion for star formation, infall or outflow of gas into galaxies, 
Lyman alpha cloud dynamics, etc. 

Most previous self-similar results are obtained as special cases. The critical values 
for a cosmological flow with = 1 and 7 = 4/3 turn out to play a special role. 

Subject headings: gravitational collapse, hydrodynamics, self-similar, shock waves 



1. Introduction 

Self-similar behavior provides an important class of unsteady solutions to the self-gravitating 
fluid equations. On one hand, many physical problems often attain self-similar limits for a wide 
range of initial conditions. On the other hand, the self-similar properties allow us to investigate 
properties of solutions in arbitrary detail, without any of the associated difficulties of numerical 
hydrodynamics. They also yield a rich class of test cases for general purpose numerical codes of 
gravitational hydrodynamics. 

In astrophysical scenarios gravity is often important, and in those cases the general form of 
similarity solutions is very restricted. In many situations the asymptotic boundary consists a cold 
flow, in which case the solution space can be qualitatively understood. We can determine the 
qualitative and asymptotic behaviour of all cold collisionless self-similar flows. 
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This paper encompasses and generalizes the calculations of previous papers, ranging from 



cosmological explosions to gravitational collapse (Ostriker and McKee 1988, Bertschinger 1985b 



Ikeuchi et al 1983, Lemos and Lynden-Bell 1989b| , etc). The following mathematical analysis yields 



more insight into the nature of solutions. 

In this paper, I will start by reviewing the self-similar ansatz. Then I identify the critical 
points, which serve as the asymptotic limit of solutions. By assuming zero pressure, we can 
qualitatively solve the problem in the velocity-density phase space. This allows us to construct 
full hydrodynamical solutions by matching a shock jump to the cold flow. After integrating the 
isothermal equations, I conclude with some suggestions for potential applications. 



2. Equations 



The goal of this paper is to solve and classify solutions to the spherically symmetric gravitating 
inviscid barytropic fluid equations: 

dp ( d d \ 1 d . 2 N /lN 

dv 1 dp Gm 

dt pdr r 2 ' 



dt 
dm 

dr 



0, (3) 
4irr 2 p. (4) 



Self-similarity allows us to reduce the self-gravitating fluid equations from partial differential 
equations (PDE) into ordinary differential equations (ODE) flSedov 1959j ). 

The general self-similar solution allows two dimensionful parameters. In the case of 
self-gravitating fluids, one is the gravitational constant with dimensions 

[G] = M~ 1 L^T~ % . (5) 

This already uniquely determines the dimensionless parametrization of any similarity solution: 

v = ~V{\), (6) 
P = Phtt(\), (7) 
p = PhCjfPW, (8) 

m = ^nr 3 p h M(X) (9) 
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where 



* 55 mx- (10> 



The second dimensionful parameter, call it A, has some dimension 

[A] = T . (11) 

We have required that it not contain the dimensions of mass, which one can eliminate using the 
gravitational constant G. Furthermore, we scaled A such that the power of length is —1. Some 
examples will be given shortly. 

We now have determined the exact form of the dimensionless independent variable A 

X = A^. (12) 
After some algebra, we obtain the equations which any self-similar solution must satisfy: 

K n ^ V n 6) +V) = 2- W (13) 

2P P' 2M 

V(V-l) + XV'(V-5) = -_-A--— (14) 

HV-S)(y-^) = 4 -2 7 -2F (15) 

Here we have integrated the mass equation, as presented by Ikeuchi, Tomisaka and Ostriker 



( Ikeuchi et al 1983 ) hereafter referred to as ITO. We thus have a third order ODE with three 
degrees of freedom for a given gas index 7 and scaling relation 5. The equations are scale invariant 
under a change A — ► a\, eliminating the freedom in the second dimensional parameter A to only 
its dimension, which determined 5. Since A is a constant of the problem, it fixes the scale if we 
express it as an integral over A. Such an example is worked out in Ostriker and Pen (1993). In 
general one can obtain this integral from the initial conditions, as I will exemplify below. 



3. Analytic Results 



We can immediately extract some information from equations (13-16|). First, we look at 
the critical points by setting all derivatives with respect to A to zero. This yields the static 
solutions which often form the boundary conditions of the general solution. From (13) we see that 
V = 2/3, unless = 0, which describes a matter dominated universe. This is an example where 
the Newtonian equations yield general relativistic exact results. Equation (O) is automatically 
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satisfied as long as P = 0. From ( |i~4| ) and (|l6|) we obtain the curious result that for 7 = 4/3 a 
solution family P = (1 — S1)S7/ 9 exists. We recover ^ = 1=^^ = as the expected limit. This 
provides a Newtonian solution for a universe filled with a gas supported by radiation pressure. 
Note from equation (||), however, that the pressure increases at fixed t as r 2 without bound, while 
decaying as t~ 4 at fixed r. The 7 = 4/3 gas seems to take a special role, as was observed by 
Bertschinger ( [Bertschinger 1985b|) for the special case of uncompensated collapse. If 7 7^ 4/3, the 
only solution is D, = 1, P = 0, corresponding to the free Hubble flow. 

We can also model a universe which is filled by a massless radiation gas, such as the early 
radiation dominated universe, by accounting for the relativistic pressure. This is achieved through 
addition of the equivalent inertial and gravitational mass of the relativistic gas. The correct 
equation of state has 7 = 1 since p = c 2 p/3. c is the speed of light, which determines A = 1/c, 
and thus 5 = 1. We require p to be constant, so P oc A~ 2 . The unit value of 5 automatically 
satisfies (|j~5|). Write the flat line element as —c 2 dt 2 + dx 2 . From the energy momentum tensor for 
a relativistic gas T^ v = (p + p/c 2 )u^u v +p/c 2 g >ll/ and the conservation law T£ u = we obtain the 
following modified continuity equation to first order in v/c: 

m' + w<>+,/*W = < 17 > 

^(p+p/c 2 )v< + A((p + p / c 2 )t ,v + j,{«) = p V3 , (is) 



We have written the 4-velocity u a = (Vc 2 — v 2 /c, v 1 ,v 2 ,v 3 )/c, and is the Newtonian potential. 
18) differs from the classical Euler equation by the appearance of p + p/c 2 in place of p. 



Then critical point now becomes V = 1/2, the general relativistic expansion law for a radiation 
dominated universe. 

Let us continue the search for critical points. At A = 00 we get another solution 
P = 0, £7 = 0, V = 0, with Q,' /Q = 2/{5\). This would correspond to an ever diminishing density 
at infinity, and might be of interest in star formation problems. The remaining choices are to let 
one of the variables diverge, with either V = which we will study in some detail, or V = —00 
which leads to a degenerate class of solutions. The limit V = 5 is especially curious, since that 
generates divergent densities moving at finite velocities. We will examine each of these possibilities 
below. 

Through a change of variables £ = ln(A), we obtain a homogeneous set of equations which do 
not contain the independent variable explicitly. We can now ask for the general form of a solution. 
One possible boundary condition at either A = or A = 00 is given by the Hubble solution. 



Furthermore, since the mass equation (16) must be positive, we have V < 5 as long as 5 > 2/3, 
otherwise V > 5. The latter case implies that the fluid speed is bounded from below. To allow 
negative velocities one needs to resort to 5 < 0. We have a discrete symmetry for A — * —A. So in 
principle one can construct solutions for t < which encounter a singularity at t = 0, similar to 
the scaling solution for texture sjSpergel and Turok 1990 . 



Now let us interpret some few choices of 5. A special choice of 4/5 corresponding to a second 
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parameter with dimensions of energy. We can express our constant A in terms of the energy E as 

A = {GE)~ 5 (19) 



which satisfies our form given in ( |Tl| ). This solution is investigated in further detail in Ostriker 
and Pen 1994 . 

We can derive the 5 = 8/9 solution as found by Bertschinger ( Bertschinger 1985b[) from a 
simple dimensional viewpoint. He starts with a top hat perturbation of mass excess Mi defined in 
his terms as 

Mi = 5 iPh rf (20) 



,2/3 

uu u 

,2/9 



at a time ij. We require GMj/r, = GMi/t/ to be a constant scale of the problem, leading us to 
the dimensional constant 



a = 1gW^ (21> 

with 5 = 8/9 as desired. We will see below how to derive the scaling more rigorously. 

An isothermal gas is characterized by its sound speed, so 5 = 1 with simply A = l/v soun( i- 
We will return to this case in more detail in section ^] If the constant A had dimensions of mass 
it would imply 5 = 2/3, which is forbidden by the mass equation ( showing that self-similar 
solutions with characteristic mass do not exist. 



4. Cold Flows 



Often we are interested in situations where the fluid is cold at large distances. Lemos and 



Lynden-Bell 1989b had looked at pressureless cold flows without shell crossing in a Lagrangian 



frame. We can now qualitatively solve all such cases in a Eulerian coordinate system. 



When we take P = 0, all solutions are limited to the $7 — V phase plane. (13- 16) become the 
simple system 



n 
n 

V(6 - 95) 
V -5 



2-3V -V 

V -5 
V(6 - 95) + 2Q 
1-V ' 



(22) 
(23) 



A dot denotes differentiation with respect to £. The homogeneous equation in two unknowns 
can be completely classified once we know the critical points. First we note that the physically 
allowed region is given by > and V(35 — 2) < 5. The critical points in this region are 
P = (Q = l, V = 2/3), P x = (O = 0, V = 0), P 2 = (O = oo, V = oo), P 3 = (0 = oo, V = 
5), P4 = (O = 0, V = 1). The last point P4 moves to P3 when 2/3 < 5 < 1. The points P3 and P4 
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differ also from the other three in the sense that the solution curve reaches those points at a finite 
value of the parameter £. 

The critical point P$ has two asymptotics leaving and two entering, dividing phase space into 
four disconnected regions. First, consider 7 7^ 4/3 and 1 > 5 > 2/3. The second condition follows 
from equation ( |i"6| ) by requiring mass to be positive. Figure [l] is characteristic of such solutions. 
It shows the £1 — V solutions near the critical point for 5 = 4/5. Two solutions originate from 
the critical Hubble flow, and two converge onto it, dividing the phase space into four regions as 
follows: 

I. Solutions are infinitely dense at the origin, with velocities going infinitely negative. One 
might try to interpret this in terms of cold matter accreting onto a black hole. Far away, density 
and velocity go to zero. These solution are very attractive to many astrophysical scenarios 
involving infall. One can see some solutions involving positive V corresponding to outflow, while 
below a certain limit all solutions involve purely negative velocities. The one remaining parameter 
freedom corresponding to the slope of the trajectory ar critical point P\ gives the modeler more 
freedom to fit self-similar evolution to astrophysical phenomena. The divergences can be prevented 
by introducing a hydrodynamic shock as discussed below. 

II. Velocities are always positive, starting at a finite value of 5, where the density is infinite, 
with the total mass M going to 0. Both V and decrease to zero at infinity. This is a collisionless 
analogy to the ITO solution of an expanding hole, except this time through a less dense medium. 
This family of solutions will also be palatable to explosive shock continuations, and may correspond 
to some real situations. 

III. The matter is constrained to a shell of finite thickness, such as the propagation of 
exploding shells through empty space (Ikeuchi et al 1983). The cumulative amount of matter m 



diverges at the outer boundary, since a finite conserved quantity is not allowed. Densities diverge 
at the inner and outer boundaries, which both move as rj, oc t s . If interpreted as expanding 
shells of matter, we are looking at a decelerating shell. So far this only holds for cold matter, 
an assumption which would not hold for explosion powered expansion. This and the next class 
may not be very physical since they require cold, divergent masses at a finite radii. Numerical 
studies show that adding a bit of pressure changes the behaviour radically. These solutions can be 
obtained by matching a shock jump, as described below. 

IV. The last region again does not seem to correspond to any common physical condition. 
The cold gas is constrained to an expanding sphere of finite radius but infinite mass, at the interior 
of which we have an accreting black hole. Again, solutions of interest may involve invoking finite 
pressure. 

Especially interesting are the asymptotes leading to and from the critical point, defining 
the border between these regions. Since the vertex point Po describes the Hubble flow, the two 
"incoming" solutions correspond to a Hubble boundary condition at spatial infinity. The left 
outgoing trajectory is the unique solution where all intensive physical quantities are bounded 
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through all of space. If the universe were created in this form, astronomers observing from the 
origin might be tricked into believing the universe to be a Hubble expanding big bang, when in 
fact densities go to zero at large distances, giving an illusion of a locally isotropic distribution. The 
right outgoing trajectory shows a cold expanding shell of matter with finite density in the center, 
diverging at some finite radius. Its existence is more a curiosity than a realistic astrophysical 
scenerio, showing the range of behaviors that similarity allows. 

We can quantify the solutions near this isolated critical point. When we expand (^) to linear 
order for small perturbations around Pi, the linear equations have eigenvalues 

h = 3^h< h = ™ 

and eigenvectors 

Vi = 

m = {"es-i )■ < 25 > 

The fact that the eigenvalues are real and have opposite sign classifies this as a hyperbolic critical 
point, which indeed divides phase space. The geometric interpretation of these vectors is the 
asymptotic tangent to the trajectories at Pi, as can be seen in figures |l|, |2]and ^. The eigenvalue 
determines the direction of the arrow. 

These vectors describe solutions connecting to the Hubble flow at infinity for the negative 
eigenvalue, and Hubble flow at the origin for the positive case. For the first case, the asymptotic 
behaviour of £1 — 1 and V — 2/3 goes as 

An application of this analysis is that Bertschinger's infall solution of 5 = 8/9 can be obtained 
by requiring Q to be monotonic and asymptotically approaching 1. The statement corresponds to 
requring a logarithmically asymptotic Hubble flow. The first component of the second eigenvector 
V2 must then be zero, leading to the correct scaling. 

Other scalings of interest include 5 = 4/5 for energy conservation, which was investigated 
in further detail by Ostriker and Pen 1994| . The isothermal case reduces to the cold flow for 




sufficiently large A, as we will show below. This makes Figure Q the correct solution in that limit. 
We note a qualitative change in the solutions, that the upper left corner has actually become a 
valid stable asymptote. 

The solution space bifurcates at 5 = 1, 2/3 and 0. The first case we consider in section p], 
and the last we see illustrated in figure ^. The bottom domain boundary is at V = 5 = —1/2, so 
this graph includes the critical point Pi. Note that the direction of the line connecting Pi to Po 
has reversed. 
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5. Shock Jump 



We now consider solutions which actually interact hydrodynamically. Typically this will occur 
when a hot gas shocks through the cold environment, or cold matter passes through an accretion 
shock onto a collapsed object. When one region of the solution (or initial condition) is very cold 
compared to the kinetic or thermal energy of another region, the previous analysis will accurately 
describe the cold domain. The hot region is matched onto it using a strong shock discontinuity. 

We could introduce a shock jump at an arbitrary point in phase space. What happens, 
however, is that such solutions either diverge, as shown for the critical point by |Ikeuchi et al| 



1983, or else collapse with very divergent dynamical values. A critical line exists in phase space, 



balanced between collapse and void, which satisfy the physical requirement of zero velocities at 



the origin. For special cases of 5 = 4/5,8/9 these solutions have previously been found (Ostriker 
and Pen 1994, and [Bertschinger 1985b| ). 



Let us examine the infall solutions. We require matter to actually decelerate and come to a 
stop as it falls in. At large distances, the matter is taken to be cold, so that the analysis of the 
previous section holds. As we follow a solution from spatial infinity along the line of decreasing A, 
we will encounter a strong accretion shock given by the Rankine Hugionot relations 

V s = 6 + ^Z±(V -S) 
7 + 1 



P, 



7-1 
2n (V -5) 2 



7 + 1 

M s = M Q . (26) 

Here a subscript of s denotes the postshock values and the preshock condition. 

This jump will move the curve up and to the left in figure ||. The shock could occur anywhere 
along the solution curve. We will break this freedom by requiring the constraint V(0) = 0. Since 
the problem is scale invariant in A, the shock jump location is not characterized by £ but by its 
location in the £1 — V plane. In order to fix the parametrization one needs to fix the definition 
of the dimensional constant in terms of an integral over the hydrodynamic quantities. After the 
shock, solutions may cross in their Q — V projection since the pressure P / 0. Figure ^ shows 
such trajectories corresponding to 5 = 4/5,7 = 5/3, connecting to the Hubble flow or a stationary 
medium respectively. We note that solutions exist not only for collapsing gas spheres, but also for 
expanding accreting ones, as the dotted line above V = represents. 

This approach reproduces the previously known solutions, and includes generalizations for 
trajectories which start with stationary media and low densities at spatial infinity. Solutions exist 
for a range in 5, corresponding to various scalings as discussed earlier. 
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If we drop the requirement of bounded velocity at the origin, e.g. in accretion problems onto 
a black hole, the shock jump location is no longer uniquely determined. A large family of solution 
exists in this case. Another way of increasing the size of the solution space is to introduce several 
shock jumps. 

The degeneracy of the solutions exists again for explosions. It relates to the fact that the self 
similar solutions remember one more parameter from the initial conditions, corresponding to the 
position of the shock jump. 



6. Isothermal Solutions 

The existence of a characteristic velocity, such as in the case of an isothermal flow, would 
lead to S = 1. In this case 7 = 1 and the pressure equation (|i~6|) is integrable as P = cqQ/X for an 
arbitrary constant cq. While the dimensionless parametrization (^) seems rather unusual, a scaling 
in terms of the sound speed e.g. v = cV yield the same equations with a change of variables 
V — > V/X. But if we nevertheless use G as the first dimensionful constant, we can integrate the 
pressure equation ( |i~5| ) to yield a new Euler equation 

v e~* (tl \ 2M 
V(V-l) + V(V-5) = -^--3y— . (27) 

In terms of the autonomous equations of ([i~6|) in the variable £, the pressure term 
exp(— £)(f2/Q — 3)/2 in (^7|) will decay exponentially as £ — > 00. For sufficiently large £, we thus 
recover the pressureless case, which we had investigated earlier. A change in behavior occurs, 
however, for the critical point V = 1, Q = 0: it is now a stable asymptotic limit, so solutions of 
this kind exist. 

An isothermal shock introduces two more free parameters: the location and strength of the 
shock, both of which are arbitrary. The strong shock conditions ( ^6|) are no longer uniquely 
determined, but we nevertheless expect the same qualitative behaviors as in the adiabatic shock, 
since the have seen the isothermal case to be identical to the general solution with 7 = 1. 



7. Conclusions 

We now have an improved understanding of the general behaviour of self-similar solutions of 
self-gravitating fluids. The solution space is of managable size, and all past solutions fit into the 
new framework. A class of new solutions has also been found under the ansatz that the flow be 
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cold at large distances. The self-similar motion for cold gas, which is equivalent to a pressureless 
fluid without shell crossing, has been solved, with some surprising asympotic behaviour in its 
— V phase space near critical points. 

Since the similarity equations were derived from the gravitating Euler equations, any 
similarity solution must be a valid physical scenario under the assumptions of a polytropic inviscid 
gas. Since we now know the possible solution space, we can compare astrophysical situations with 
the given templates and decide on the applicability of the self-similar ansatz. 
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Figure Captions 

Figure ||: Cold flow solution for S = 4/5. Arrows indicate direction of increasing £. The 
heavy lines describe the four asymptotic flows with a Hubble boundary condition. 

Figure ^: Cold flow solution for 5 = —1/2. Here the velocity field is bounded from below 
because 5 < 2/3. 

Figyure |3|: Infall solutions for 5 = 4/5 The heavy line marks the boundary in phase space 
where a shock jump connects onto the cold flow. The solid lines represent a solution to the Hubble 
flow at infinity, while the dashed line shows a solution from a stationary medium. The dotted line 
is an expanding collapse solution. 

Figure |j: Isothermal flow solution with 5=1. Any isothermal solution must converge to 
these trajectories for sufficiently large A. 
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5 = 4/5 




0.5 1 1.5 2 

Q 

Fig. 1. — Cold flow solution for S = 4/5 




Fig. 3. — Infall solutions for (5 = 4/5 
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6=1 




0.5 1 1.5 

Q 

Fig. 4. — Isothermal flow solution with 5 = 1 

Any isothermal solution must converge to these trajectories for sufficiently large A. 



